Application of the inverse Hamiltonian method to Hartree-Fock-Bogoliubov 

calculations 
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We solve the Hartree-Fock-Bogoliubov (HFB) equations for a spherical mean field and a pairing 
potential with the inverse Hamiltonian method, which we have developed for the solution of the 
Dirac equation. This method is based on the variational principle for the inverse Hamiltonian, and 
is applicable to Hamiltonians that are bound neither from above nor below. We demonstrate that 
the method works well not only for the Dirac but also for the HFB equations. 
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Pairing correlations between nucleons play a very im- 
portant role in open shell nuclei [l|, [^. Hartree-Fock- 
Bogoliubov (HFB) theory is a powerful method which 
treats these correlations in a self-consistent way in the 
framework of a single generalized Slater determinant of 
independent quasi-particles [2|-|j|. The method has been 
widely used in recent years for the study of the structure 
of neutron rich nuclei far from stability up to the neu- 
tron drip line, where the coupling to the continuum has 
a great influence. 

The HFB-equations are a set of coupled differential 
equations. In many cases they have been solved by an 
expansion in terms of a finite set of basis functions, as 
for instance the eigenfunctions of a harmonic oscillator 
or a Woods-Saxon potential. Although this method has 
been used successfully for many investigations in the lit- 
erature it has its limitations, (i) The convergence with 
the number of basis functions depends on the parame- 
ters of the basis and the optimization of this parameters 
is often very complicated, (ii) in many cases, in particu- 
lar for two-dimensional (2D) and three-dimensional (3D) 
calculations in heavy nuclei, the dimension of the ma- 
trices becomes extremely large. Since the CPU time for 
one diagonalization grows with the cube of this dimen- 
sion this is connected with considerable numerical efforts, 
(iii) In each step of the iteration the corresponding ma- 
trix elements of the two-body interaction in terms of the 
basis functions requires in addition much CPU time, (iv) 
the treatment of the continuum is connected with specific 
difficulties in particular in the case of neutron halos. In 
summary these methods do not exploit fully the advan- 
tages of powerful zero-range interactions as for instance 
the Skyrme energy density functional. These advantages 
can only be exploited fully in coordinate space. For this 
reason the Orsay group has introduced already in the 
eighties for the solution of the Hartree Fock (HF) equa- 
tions the imaginary time method [5|, |6| . Starting with an 
initial HF wave functions Ivt^^)) the solution is obtained 
in iterative steps with infinitesimal step size t as 
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HF equation. An important fact is that the spectrum of 
the HFB equation is bound neither from above nor from 
below (see Fig. [TJa)) because of the coupling between 
particle creation and annihilation parts in the quasi- 
particle operators of the Bogoliubov transformation 0] . 
This inhibits a direct application of the imaginary time 
method to HFB, which has been successfully employed 
in self-consistent mean field calculations in the coordi- 
nate space representation [g, |7[- That is, if the imaginary 
time evolution is naively applied, the iterative solution in- 
evitably dives into the quasi-particle negative continuum. 
To avoid this problem the rather coinplicated two-basis 
method has been introduced in Refs. [3,|8| where, in each 
step of the iteration, the HFB equations are solved by ex- 
pansion of the quasi-particle wave functions in a HF-basis 
calculated by the imaginary time method on a 3D mesh 
in coordinate-space. In Ref. Q the HFB equation on a 
3D mesh has been solved in the canonical basis. 

It should be mentioned that the diagonalization of 
the huge matrices in the basis expansion method can be 
avoided by the gradient method introduced in Ref. [ifll, 
which is also applicable to the solution of the HFB equa- 
tion with a spectrum not bound from below. Here the 
wave function in the next step of the iteration is expressed 
in terms of the Thouless theorem [21 
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Obviously this method is limited to Hamiltonians H with 
a spectrum bounded from below, as the non-relativistic 



where |v[/(")) ig the vacuum with respect to the quasi- 
particle operators a\. For infinitesimal r only 2-quasi- 
particle states with positive energy Ek + E^i are admixed 
in the next step of the iteration. This method has been 
successfully used in the literature for the solution of the 
HFB equations in an oscillator basis [11[ and it has been 
also applied for a variation after projection (VAP) in 
Refs. [i3:ll3l- Of course, the matrix elements H^^, have 
to be calculated here in the corresponding quasi-particle 
basis. Therefore the locality of the Skyrme interaction 
and the quasi-locality of the kinetic energy in coordinate 
space cannot be exploited in this method. 

We therefore propose in this report for the solution of 
the HFB equations a different method keeping in mind 



that one is confronted with the same problem in the so- 
lution of the Dirac equation because the corresponding 
spectrum has the Dirac sea states down to the negative 
infinity as well as the positive energy states up to the pos- 
itive infinity. This leads to a breakdown of variational 
calculations which has long been known in the field of 
relativistic quantum chemistry under the name of "varia- 
tional collapse" , and there has been a number of prescrip- 
tions proposed to avoid it [l4l - [l8| . Recently, newly devel- 
oped methods for iterative solutions of a Dirac equation 
are introduced by Zhang et al. [l^ [20| and by Hagino 
and Tanimura |2l| in the nuclear physics context. In 
Ref . [21| , based on the idea of Hill and Krauthauser [13] 
a novel method has been developed, which is called "in- 
verse Hamiltonian method" , for relativistic mean field 
calculations in the coordinate space representation. In 
this method a variational principle is applied to the Her- 
mitian operator l/(ilDirac — W) instead of the Hamilto- 
nian -ffoirac itsclf. Here, VT is a real number which is 
set between the Fermi sea and the Dirac sea. In con- 
trast to some other methods [H, [la, [13, 1201 , it is rela- 
tively straightforward to apply our method not only to 
the Dirac equation but also to other eigenvalue problems 
with unbound operators, such as HFB. In this paper, we 
apply the inverse Hamiltonian method to a HFB equa- 
tion in the coordinate space representation and show that 
the equation can be solved successfully without the vari- 
ational collapse. 

In HFB calculations one usually needs to obtain several 
lowest positive energy quasi-particle states only. That 
is, one only needs the states ^1,^2, ■■■ associated with 
eigenvalues Ei,E2,... (see Fig. 1(a)) As is seen in Fig. 
[TJb), these states come to the top of the spectrum of 
1/{H — W) if W is set between the positive and negative 
spectra. A variational principle for 1/{H — W) is that a 
maximization oi {{H — W)~^) leads to the desired quasi- 
particle state wave functions ^^. Our method maxi- 
mizes {{H — W)~^) based on the relation [21j 
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where \ip^'^') is an arbitrary wave function which is not 
an eigenfunction of the Hamiltonian H and VF is a real 
constant between Ei and E^i. All the states below W 
damp out and only i/)! which is just above W survives in 
the limit n — >■ 00. 

In practice, the wave function ■0("+i) is evolved for a 
small step AT from -0^"^ as 
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To this end, we have to solve a large sparse linear equa- 
tion 
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in order to invert the Hamiltonian. We here employ an 
iterative method for linear systems, that is, the conjugate 



(a) 



.. . E 2 E ^ t J E2 ■ . . 
H OO O O I • • ••!- 



->E 



(b) 



O O OO h 



1^ 

V(-l-W) 



l/(Ei-W) 

-I • • • i > 






FIG. 1: (Color Online) Spectra of (a) a quasi-particle Hamil- 
tonian H itself and (b) the inverse of the Hamiltonian 1/{H — 
W). A is the chemical potential. The bound states of pos- 
itive and negative energies are indicated by solid and open 
circles, respectively. The continuum states are represented by 
the thick solid lines. The energy shift W is taken between the 
positive and negative spectra. The eigenvalues are labeled by 
an integer k such that E^k = —Ek- 



gradient normal residual (CGNR) method 22]. This is 
one of the Krylov subspace methods for sparse linear sys- 
tems [22, [231 • CGNR solves a linear system Ax = 6 by 
applying the conjugate gradient method to an equivalent 
system A'' Ax = A^b. 

When the mean field and pairing potentials are local 
and spherical, a quasi-particle wave function is given by 
the form 
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Here, %,jm is a spherical spinor defined by ^j„i = 
I]m,m'(^"^5™'b'"^)^^'«Xm' : whcrc Y^^ and Xm' are 
spherical harmonics and spin wave function, respectively. 
The HFB equation in the coordinate space then reduces 
to a the radial equation 
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where h is the mean field Hamiltonian, A is the chemical 
potential, and A(r) is the pairing potential. Following 
Refs. [23jl23) we use a phenomenological Woods-Saxon 
type potentials, which simulates medium-heavy neutron- 
rich nuclei around ^^Ni, for the mean field and the pairing 
potentials. The potential v{r) in the mean field Hamil- 
tonian h and the paring potential A(r) are thus taken 
as 
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with vo = 
fm, and a 



38.5 MeV^^= = 14 MeV-fm^, Rq = 5.63 
0.66 fm 2J: l25|- The strength of pairing 



potential Aq is determined so that the average pairing 
gap A defined by ^2d\ 
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is equal to 1.0 MeV. The chemical potential A is fixed to 
A = —0.5 MeV in the present calculation. We solve Eq. 
([7|) by discretizing the radial coordinate r with mesh size 
Ar, and imposing the box boundary condition. The sec- 
ond derivative of ip at the ith mesh point is approximated 
by 3-point difference: ijj'l = {ipi+i — 2^pi + ?/)i_i)/(Ar)^. 
Let us now apply the inverse Hamiltonian method 
and numerically solve the HFB equation, Eq. ([7]). We 
also solve the equation exactly by directly diagonaliz- 
ing the coordinate space Hamiltonian. by the Runge- 
Kutta method. The parameters of the inverse Hamil- 
tonian method are set W = 0.1 MeV and AT = 10 
MeV. The excited states are also calculated simultane- 
ously by orthogonalizing a set of wave functions at every 
step of iteration. The radial coordinate is discretized up 
to Tmax = 30 fm with Ar = 0.1 fm. Initial quasi-particle 
wave functions are taken to be a Gaussian form 
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where £ is the orbital angular momentum and Nk is an 
appropriate normalizing constant. The width parameter 
of the Gaussian bk is taken as bk = 2.5 x 1.05*^"^ fm, 
(fc = l,2,...). 

Let us first discuss the convergence properties of the 
energy (H) and the expectation value of the inverse of 
Hamiltonian {{H — W)~^) for the lowest si/2 quasi- 
particle state. In Fig. [2l we show the evolution of the 
two quantities as functions of the number of iteration 
steps. As is observed in Ref. [2l| for a Dirac equation, 
{{H — W)~^) converges monotonically up to a certain 
value as the iteration step increases. At the same time, 
(H) converges to the lowest Si/2 eigenvalue, E — 0.424 
MeV. 

In Table HI we show quasi-particle energies and occu- 
pation probabilities v'^ for the three lowest Si/2 states in 
comparison with the exact values which are obtained by 
diagonalizing the Hamiltonian. The occupation probabil- 
ities are defined in terms of quasi-particle wave function 
by 
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The agreement is perfect both in the energies and the oc- 
cupation probabilities for the digits shown in the table. 
Fig. [3] shows comparisons of wave functions of the three 
Si/2 states. The dashed lines show the exact wave func- 
tions, whereas the solid lines show the wave functions 
obtained with the inverse Hamiltonian method. The 
left and right panels show the upper component Uk{r) 
and the lower component Vk (r) of a quasi-particle wave 
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FIG. 2: (Color Online) Covergence properties of (a) the en- 
ergy expectation value (H) and (b) the expectation value of 
the inverse of Hamiltonian {{H — W)~^) for the lowest Si/2 
quasi-particle state. The energy shift and the step size of T 
are taken to be W = 0.1 MeV and AT = 10 MeV, respec- 
tively. 



TABLE I: A comparison between the exact calculations and 
the inverse Hamiltonian method for the three lowest Si/2 
quasi-particle energies E and occupation probabilities v^. 
The exact values are calculated by diagonalizing the real space 
Hamiltonian. 



E (MeV) 




vl 


exact inv. H method 


exact 


inv. H method 


0.42414 0.42414 


0.5574 


0.5574 


1.0383 1.0383 


3.972 X 10" 


'2 3.972 X 10"^ 


2.3063 2.3063 


9.689 X 10" 


-•' 9.689 X 10-^ 



function, respectively. As is seen in Fig [3J the inverse 
Hamiltonian method reproduces the wave functions al- 
most identically to the exact ones for both the bound 
state and the excited continuum states. We have also ob- 
tained the other s-wave states with an accuracy as high 
as the lower states shown in Table U and Fig. [3l 

We have checked the performance of the inverse Hamil- 
tonian method for other angular momentum quantum 
numbers and confirmed that the method solves the HFB 
equation as accurately as for the Si/2 states. It is ap- 
parent that the inverse Hamiltonian method gives prac- 
tically the exact solutions of the HFB equation in the 
coordinate space representation and is safe against the 
variational collapse. 

In summary, we have discussed the numerical perfor- 
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FIG. 3: (Color Online) Comparisons of wave functions for 
the three lowest Si/2 states. The left and right panels show 
the upper {Uk{r)) and the lower (Vfe(r)) components of quasi- 
particle wave function, respectively. The exact wave functions 
are shown with the dashed lines and the ones obtained by the 
inverse Hamiltonian method are drawn with the solid lines. 



mance of the inverse Hamiltonian method for a HFB cal- 
culation. While the method has been developed for solv- 
ing Dirac equations, we have shown that it can almost 
exactly solve a coordinate space HFB equation as well 
with spherical mean field and pairing potentials without 
variational collapse. An obvious future work is the appli- 
cation of the method to self-consistent HFB calculations 
on 3D mesh. A work in this direction is now in progress. 
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